% MODEL.m
%
% Sets up the Stationary version of Ireland's (2007) New Keynesian model that we use for the applications in our paper. 
%
%                        1       2     3     4         5              6          7       8        9       10 
% Variables in Y(t) = [ yhat(t) pi(t) r(t) ghat(t) E_t(yhat(t+1)) E_t(pi(t+1)) ahat(t) ehat(t) ln(z(t)) zhat(t)]'
% 
%
%                           1      2    3     4
% Variables in eps(t) = [ ea(t) ee(t) ez(t) er(t) ]'
%
%
% Variables in eta(t) = [etay(t) etapi(t)]'
%
% Date Created: 7-April-2008
% Last Modified: 1-May-2008


%% Load initial and final parameters.
param

%% Checks

if or(ta < 2, up_to < T+ta) 
    disp('Possible Error: increase time of annoucement, ta')
    disp('Possible Error: increase time of simulation, up_to')
    return
end

%% Load Dimensions
nvars = 10;
nshocks = 4;
nfcast = 2;

%% Load Initial Matrices
% Note: row indicates equation number and column indicates variables number
GAM0 = zeros(nvars);
GAM0(1,1) = 1; 
GAM0(1,3) = 1/sigma;
GAM0(1,5) = -1;
GAM0(1,6) = -1/sigma; 
GAM0(1,7) = -(1-rho_a)/sigma;
GAM0(2,1) = -psi*sigma;
GAM0(2,2) = 1+beta*alpha;
GAM0(2,6) = -beta;
GAM0(2,8) = 1;
GAM0(2,10) = psi;
GAM0(3,1) = -rho_y;
GAM0(3,2) = -rho_pi;
GAM0(3,3) = 1;
GAM0(3,4) = -rho_g;
GAM0(4,7) = 1;
GAM0(5,8) = 1;
GAM0(6,9) = 1;
GAM0(7,1) = -1;
GAM0(7,4) = 1;
GAM0(8,9) = -1;
GAM0(8,10) = 1;
GAM0(9,1) = 1;
GAM0(10,2) = 1;

GAM1 = zeros(nvars);
GAM1(2,2) = alpha;
GAM1(3,3) = rho_r;
GAM1(4,7) = rho_a;
GAM1(5,8) = rho_e;
GAM1(6,9) = rho_z;
GAM1(7,1) = -1;
GAM1(9,5) = 1;
GAM1(10,6) = 1;

C = zeros(nvars,1);
C(1,1) = -(1/sigma)*log(beta);
C(2,1) = (1+beta*alpha-alpha-beta)*pitarget;
C(3,1) = (1-rho_r)*(pitarget-log(beta)) - rho_pi*pitarget;
C(6,1) = (1-rho_z)*log(z);
C(8,1) = -log(z);

PSI = zeros(nvars,nshocks);
PSI(3,4) = 1;
PSI(4,1) = 1;
PSI(5,2) = 1;
PSI(6,3) = 1;

PPI = zeros(nvars,nfcast);
PPI(9,1) = 1;
PPI(10,2) = 1;

%% Load Final Matrices
GAM0n = zeros(nvars);
GAM0n(1,1) = 1; 
GAM0n(1,3) = 1/sigma_n;
GAM0n(1,5) = -1;
GAM0n(1,6) = -1/sigma_n; 
GAM0n(1,7) = -(1-rho_a_n)/sigma_n;
GAM0n(2,1) = -psi_n*sigma_n;
GAM0n(2,2) = 1+beta_n*alpha_n;
GAM0n(2,6) = -beta_n;
GAM0n(2,8) = 1;
GAM0n(2,10) = psi_n;
GAM0n(3,1) = -rho_y_n;
GAM0n(3,2) = -rho_pi_n;
GAM0n(3,3) = 1;
GAM0n(3,4) = -rho_g_n;
GAM0n(4,7) = 1;
GAM0n(5,8) = 1;
GAM0n(6,9) = 1;
GAM0n(7,1) = -1;
GAM0n(7,4) = 1;
GAM0n(8,9) = -1;
GAM0n(8,10) = 1;
GAM0n(9,1) = 1;
GAM0n(10,2) = 1;

GAM1n = zeros(nvars);
GAM1n(2,2) = alpha_n;
GAM1n(3,3) = rho_r_n;
GAM1n(4,7) = rho_a_n;
GAM1n(5,8) = rho_e_n;
GAM1n(6,9) = rho_z_n;
GAM1n(7,1) = -1;
GAM1n(9,5) = 1;
GAM1n(10,6) = 1;

Cn = zeros(nvars,1);
Cn(1,1) = -(1/sigma_n)*log(beta_n);
Cn(2,1) = (1+beta_n*alpha_n-alpha_n-beta_n)*pitarget_n;
Cn(3,1) = (1-rho_r_n)*(pitarget_n-log(beta_n)) - rho_pi_n*pitarget_n;
Cn(6,1) = (1-rho_z_n)*log(z_n);
Cn(8,1) = -log(z_n);

PSIn = zeros(nvars,nshocks);
PSIn(3,4) = 1;
PSIn(4,1) = 1;
PSIn(5,2) = 1;
PSIn(6,3) = 1;

PPIn = zeros(nvars,nfcast);
PPIn(9,1) = 1;
PPIn(10,2) = 1;

[n l] = size(PSI);
%% Intial condition
% Start from the steady state of the original system
[S0_old, S1_old, S2_old, S3_old, unique_old] = Smats(C, GAM0, GAM1, PSI, PPI);
y0 = inv(eye(n) - S1_old)*S0_old;

%% Final steady state
[S0_new, S1_new, S2_new, S3_new, unique_new] = Smats(Cn, GAM0n, GAM1n, PSIn, PPIn);
y0_new = inv(eye(n) - S1_new)*S0_new;

%% FiniteSmats
n1 = 6;
n2 = 2;
nleads = 1;

nperiods = T;
eps = zeros(l,nperiods); %
[yseq, S0, S1, S2, S3, unique] = finitere_smats(n1,n2,nleads,nperiods,GAM0n,GAM1n,Cn,PSIn,PPIn,GAM0,GAM1,C,PSI,PPI,y0,eps);

%% Compute longer simulation to check convergence to new steady state

time = 1:up_to;
ysim = zeros(n,up_to);

for t=1:(ta-1)
    ysim(:,t) = y0;
end
% ysim(:,1) = y0;
ysim(:,ta:T+(ta-1)) = yseq;

for t = T+ta:up_to
    ysim(:,t) = S0 + S1* ysim(:,t-1);
end

%% Plot the response
output_resp_d = ysim(1,:);
g_resp_d = ysim(1,2:end)-ysim(1,1:end-1); 
g_resp_d = [0, g_resp_d];

inflation_resp_d = ysim(2,:);
r_resp_d = ysim(3,:);
inflationE_resp_d = ysim(6,:);
realr_resp_d = ysim(3,:) - ysim(6,:);

% In annualised per cent
inflation_resp_d = ((exp(inflation_resp_d) - 1))*400;
r_resp_d = ((r_resp_d + 1).^4 - 1)*100;
inflationE_resp_d = ((exp(inflationE_resp_d) - 1))*400;
realr_resp_d = (exp(realr_resp_d) - 1)*400;
g_resp_d = ((g_resp_d + 1).^4 - 1)*100;

if doplots_stationary
subplot(3,2,1), plot(time,output_resp_d,'b','LineWidth',lw),title('Output'),grid,set(gca,'xtick',[ta, T+ta])
subplot(3,2,2), plot(time,g_resp_d,'b','LineWidth',lw),title('Growth Rate of Output'),grid,set(gca,'xtick',[ta, T+ta])
subplot(3,2,3), plot(time,inflation_resp_d,'b','LineWidth',lw),title('Inflation'),grid,set(gca,'xtick',[ta, T+ta])
subplot(3,2,4), plot(time,inflationE_resp_d,'b','LineWidth',lw),title('Inflation Expectations'),grid,set(gca,'xtick',[ta, T+ta])
subplot(3,2,5), plot(time,r_resp_d,'b','LineWidth',lw),title('Nominal Interest Rate'),grid,set(gca,'xtick',[ta, T+ta])
subplot(3,2,6), plot(time,realr_resp_d,'b','LineWidth',lw),title('Real Interest Rate'),grid,set(gca,'xtick',[ta, T+ta])
end
